We need the sf package:
> if (! "sf" %in% rownames(installed.packages())) install.packages("sf")
> library(sf)
Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
The coordinates of the houses of the Hanam cohort are in a zipping shapefile that can be downloaded as so:
> download.file("https://www.dropbox.com/s/yyg6gx0nycpv60b/HH%20point%20Hanam.zip?raw=1", "hanam.zip")
Once downloaded, unzip the file:
> unzip("hanam.zip")
Once unzipped read the shapefile:
> hanam <- sf::st_read("Household_point.shp")[, "Name"]
Reading layer `Household_point' from data source `/Users/choisy/Dropbox/aaa/projects/Hanam/Household_point.shp' using driver `ESRI Shapefile'
Simple feature collection with 263 features and 22 fields
geometry type: POINT
dimension: XY
bbox: xmin: 105.9174 ymin: 20.48691 xmax: 105.9356 ymax: 20.51665
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
And clean the disk
> file.remove(grep("^Household", dir(), value = TRUE))
[1] TRUE TRUE TRUE TRUE TRUE TRUE
and
> file.remove("hanam.zip")
[1] TRUE
Here the data, which is of sf class
> hanam
Simple feature collection with 263 features and 1 field
geometry type: POINT
dimension: XY
bbox: xmin: 105.9174 ymin: 20.48691 xmax: 105.9356 ymax: 20.51665
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
First 10 features:
Name geometry
1 H001 POINT (105.9274 20.49189)
2 H002 POINT (105.9286 20.49196)
3 H003 POINT (105.9314 20.49663)
4 H004 POINT (105.9274 20.49191)
5 H005 POINT (105.9303 20.49529)
6 H006 POINT (105.9329 20.4953)
7 H007 POINT (105.9296 20.49722)
8 H009 POINT (105.9266 20.49759)
9 H010 POINT (105.9296 20.49672)
10 H011 POINT (105.9306 20.49684)
The coordinates of the houses can be extracted from that object and converted into a data frame as so:
> houses <- cbind(house = hanam$Name, as.data.frame(st_coordinates(hanam)))
which gives:
> head(houses)
house X Y
1 H001 105.9274 20.49189
2 H002 105.9286 20.49196
3 H003 105.9314 20.49663
4 H004 105.9274 20.49191
5 H005 105.9303 20.49529
6 H006 105.9329 20.49530
We need the OpenStreetMap package:
> if (! "sf" %in% rownames(installed.packages())) install.packages("OpenStreetMap")
> library(OpenStreetMap)
Registered S3 methods overwritten by 'ggplot2':
method from
[.quosures rlang
c.quosures rlang
print.quosures rlang
Below is the code you’d need to download the tiles from 3 different types of maps:
> upperleft <- c(20.525681, 105.905383)
> lowerright <- c(20.463343, 105.969019)
> bing <- openmap(upperleft, lowerright, type = "bing", minNumTiles = 20)
> osm <- openmap(upperleft, lowerright, type = "osm", minNumTiles = 20)
> esri <- openmap(upperleft, lowerright, type = "esri", minNumTiles = 20)
But the downloading would be twice as fast if you downlaod the objects directly from here:
> download.file("https://www.dropbox.com/s/26y2pgouodj4rpe/bing.rds?raw=1", "bing_file.rds")
> download.file("https://www.dropbox.com/s/rd34k3ixwthzg9g/esri.rds?raw=1", "esri_file.rds")
> download.file("https://www.dropbox.com/s/827ze3xr0orbab0/osm.rds?raw=1" , "osm_file.rds")
> bing <- readRDS("BING.rds")
> esri <- readRDS("ESRI.rds")
> osm <- readRDS("OSM.rds")
> for(file in paste0(c("bing", "esri", "osm"), "_file.rds")) file.remove(file)
The function that makes the map:
> plot_hh <- function(map, points) {
+ require(sf)
+ plot(st_geometry(st_transform(points, map$tiles[[1]]$projection)))
+ plot(map, add = TRUE)
+ plot(st_geometry(st_transform(points, map$tiles[[1]]$projection)), add = TRUE, col = "red")
+ }
OSM map:
> plot_hh(osm, hanam)
ESRI map:
> plot_hh(esri, hanam)
BING map:
> plot_hh(bing, hanam)
You can calcultate distances between 2 points with the sf::st_distance(). For that, you need to install the lwgeom package:
> if (! "lwgeom" %in% rownames(installed.packages())) install.packages("lwgeom")
Once this is done, the calculation of the distance between the first 2 houses is done as so:
> sf::st_distance(hanam[1, ], hanam[2, ])
Units: [m]
[,1]
[1,] 129.5865
You can get the distance matrix directly like so:
> distmat <- sf::st_distance(hanam)
Clustering:
> hc <- hclust(as.dist(distmat))
Then, you can generate clusters, depending on a distance cut-off. For example with a cut-off of 1000 m, your clusters would be:
> cutree(hc, h = 1000)
[1] 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 2 1 1 1 2 1 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 2 1 2 1 1 1 1
[67] 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 4 4 4 4 4 4 4 4 4 4 4 4 2 1 4 4 4 4 4 4 4 5 3 3 6 6 3 3 3 3 3 3 6 3
[133] 6 6 3 5 3 6 6 3 6 3 6 6 5 5 6 6 6 3 6 3 6 3 3 6 6 3 6 5 3 3 3 3 3 6 3 3 3 3 3 6 3 3 6 6 3 3 3 6 6 6 6 6 6 6 6 6 6 7 7 8 8 8 8 8 8 8
[199] 3 8 8 7 8 7 8 5 8 8 8 3 8 8 8 8 8 8 7 7 7 7 7 7 5 8 7 8 9 9 9 1 1 1 3 9 6 6 6 9 6 3 1 5 3 6 3 9 3 6 9 9 9 9 9 5 1 5 1 2 1 1 6 6 6
You can merge this information with the house identifier in a data frame as so:
> clusters <- data.frame(house = hanam$Name, cluster = cutree(hc, h = 1000))
Which gives:
> head(clusters)
house cluster
1 H001 1
2 H002 1
3 H003 1
4 H004 1
5 H005 1
6 H006 2